% 1) construct m = {m(tj)}, a matrix where the rows are funds and the
% columns are the SDF corresponding to each of the J cash flows of the fund
% (note: these are padded with zeros to make a balanced matrix)
% 2) compute gt = m*r', where r={r(tj)} i.e. matrix with fund cash flows
% 3) compute dg = 

% r = t x nq x J
% x = t x nx x J 

function [gt, dg] = fundMCexpaff(b,x,r,p,h,fast)   %h is only there to be consistent with other functions that need h

t = size(r,1);      % # funds
nq = size(r,2);     % number of observed assets to be priced (e.g., fund mkt rf) 
nx = size(x,2);     % # pricing factors (e.g. const mkt) 

J = size(x,3);      % max # cash flows in fund 

mm = zeros(t,1,J);

x(isinf(x)) = 0;    % these are the periods post end-of-fund that were padded with zeros. 
x(isnan(x)) = 0;    % ... and VFfund_CF always zero for those obs


for j = 1:J         % # cycle through all fund cash flow periods   
    mm(:,1,j) = exp( squeeze(x(:,:,j)) * b );      % sdf t x 1 x J
end
m = repmat(mm,[1 nq 1]); %to have same dim as cash flow array r: t x nq x J
gt = sum(bsxfun(@times,r,m),3) - p;    % t x nq

 
dg = zeros(nx,nq); 

if fast == 0 || nargin < 6
   dgt = zeros(nx,nq);  
   for i = 1:t 
      if nx > 1  
         adddg = squeeze(x(i,:,:))*(squeeze(m(i,:,:)).*squeeze(r(i,:,:)))';
      else 
         adddg = permute(x(i,:,:),[1 3 2])*(squeeze(m(i,:,:)).*squeeze(r(i,:,:)))'; 
      end
      dgt = dgt+adddg;
   end
   
   dg = dgt/t;

end

%just for debugging
%global sdf
%sdf = squeeze(mm); 
%global gs
%gs = gt; 

%{
size(dg)
size(cov(gt))
W = eye(3);
inv(dg*dg')
e = eig(dg*W*dg');
e(1)
e(2) 
dg*W*cov(gt)*W*dg'
W(3,3) = 0.00001;
dg*W*dg'
e= eig(dg*W*dg');
e(1)
e(2) 
ddg = dg(:,1:2); 
e= eig(ddg*ddg');
S = cov(gt);
e(1)
e(2)
(1/t)*inv(ddg*ddg')*ddg*S(1:2,1:2)*ddg'*inv(ddg*ddg') 
%}
